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We present several improvements of the infinite matrix product state (iMPS) algorithm for finding 
ground states of one-dimensional quantum systems with long-range interactions. As a main new 
ingredient we introduce the superposed multi-optimization (SMO) method, which allows an efficient 
optimization of exponentially many MPS of different length at different sites all in one step. Hereby 
the algorithm becomes protected against position dependent effects as caused by spontaneously 
broken translational invariance. So far, these have been a major obstacle to convergence for the 
iMPS algorithm if no prior knowledge of the systems translational symmetry was accessible. Further, 
we investigate some more general methods to speed up calculations and improve convergence, which 
might be partially interesting in a much broader context, too. As a more special problem, we also 
look into translational invariant states close to an invariance braking phase transition and show how 
to avoid convergence into wrong local minima for such systems. Finally, we apply the new methods 
to polar bosons with long-range interactions. We calculate several detailed Devil's Staircases with 
the corresponding phase diagrams and investigate some supersolid properties. 
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I. INTRODUCTION 

A numerical method for the simulation of large quan- 
tum systems needs to meet two requirements: i) an 
ansatz suitable for the problem in question, and ii) ef- 
ficient algorithms to find the (at least nearly) optimal 
solution within the chosen ansatz. For one-dimensional 
quantum systems on lattices, the currently most powerful 
numerical tools are matrix product states (MPS) based 
algorithms including the density matrix renormalization 
group (DMRG) [THE]. Their primary limitation is given 
by the amount of entanglement they can handle. Several 
extensions of MPS have been conceived to overcome this 
restriction as, (e.g. [SHE]), but for most practical appli- 
cations MPS are still the first choice. This is mainly due 
to the performance of the underlying optimization rou- 
tines. Although the general task to find ground states is 
known to be NP hard [9, 10_, commonly used algorithms 
seem to have no problem to attain optimal MPS solutions 
within computer precision for a plenitude of physical rel- 
evant systems. Nevertheless, for some physical systems 
of interest these algorithms still face severe difficulties. 

In this paper, we treat such problematic cases given 
by ground states of infinite spin chains with long-range 
interactions. An increasing interest in reliable numerical 
methods for these states is, e.g., triggered by the excel- 
lent experimental control of ultracold gases and the pos- 
sibility to realize systems with long-range dipole-dipole 
interactions such as Rydberg atoms or polar molecules 
|llfH5] . Although these systems are of finite size, one is 
often interested in the thermodynamical limit (i.e., infi- 
nite systems) for a better insight. 

Different strategies are known for the numerical study 
of ground states in the thermodynamic limit. One might 
try to extrapolate results from a series of increasingly 
large finite systems [rH - TTB] or directly construct the in- 
finite state itself. The latter is, e.g., done by the infi- 
nite time-evolving block decimation (iTEBD) algorithm 



|17j . which is based on an explicit translational invari- 
ant ansatz. This ansatz is quite elegant for interactions 
which are restricted to nearest neighbors, while it gets 
impractical for long-range interactions. 

A comfortable way to incorporate long-range interac- 
tions is to encode them into a matrix product operator 
(MPO) [18-21 . This concept can be integrated in an 
infinite matrix product state (iMPS) algorithm |22] I23|. 

The basic idea of this iMPS algorithm is to obtain the 
ground state of an infinite system as the fixed point of 
a constantly growing finite state by inserting iteratively 
new sites into its middle until convergence is reached. A 
major disadvantage of this approach is that the algorithm 
generally fails to converge if the ground state has a non- 
trivial translational symmetry. 

In this paper, various extensions to the basic iMPS 
algorithm are presented. As a central element, the su- 
perposed multi-optimization (SMO) is introduced (sec- 
tion IIIB), which provides a remedy for the just men- 
tioned convergence problem. Here, the key idea is to 
join the optimization of exponentially many MPS in a 
superposition and solve it efficiently. Due to this super- 
posed optimization, the effective overall problem becomes 
translational invariant again and poses no longer a hin- 
drance for convergence. 

We also introduce several improvements which work 
independently of the SMO method. As such, we present 
two different modifications of the MPS optimization rou- 
tine: On the one hand, we use physical insight for systems 
close to translational invariance breaking phase transi- 
tion to suggest a method to reduce the danger of be- 
ing trapped into a local minimum of the energy (see- 
On the other hand, we provide a more tech- 
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nical discussion how to recycle information from previous 
optimizations to speed up calculations (section IV C ) . As 
a part of these considerations we provide a simple im- 
plementation of a Davidson like algorithm (24] based on 
information from previous optimizations (appendix |D|). 
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This implementation is not bonded to the MPS frame- 
work and hence might be used in much broader context. 

All algorithms are described in depth, such that read- 
ers who are willing to reproduce our results should find 
all the information needed for successful programming. 
As a consequence, readers who just like to understand 
the crucial ideas might find the amount of algorithmic 
details far beyond their interest. Having these two types 
of readers in mind, we partitioned the material accord- 
ing to its level of detail into different chapters, such that 
entire passages can be omitted. 

In section [TTJ we review the basic concepts and the 
iMPS algorithm as presented in [22 . Readers who feel 
save to skip this part find in Fig. [T] and Fig. [2] a picto- 
rial description of all symbols used in this section. In 
Section |III| the main new concepts of this paper are pre- 
sented, above all the SMO method. Changes of the al- 
gorithm are kept to a minimum in here in contrary to 



section IV where several algorithmic improvements and 
their numerical realization are presented in detail. Each 
subsections of section |IV| contains an individual topic, 
which is presented in the first lines. These subsections 
can be skipped without danger of losing the ability to un- 
derstand the rest of the paper. A slight exception might 
be section |IV C[ in which the concept of iterative eigen- 
vector solver based on subspace projections is reviewed. 
Familiarity with this concept is assumed in section |IV D| 
|IV E| and appendix |B 1| In section [V] we apply our new 
algorithm to a system of polar bosons with long-range in- 
teractions. Detailed calculations of Devil's Staircases and 
phase diagrams are shown and a supersolid like phase 
is investigated. Finally, the paper is complemented by 
an appendix, into which several details have been out- 
sourced. 



tion is needed. A matrix product state (MPS) [TJ [3J H] 
consists in the ansatz 
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where we used the Einstein summation convention. For 
a general exact quantum state exponentially growing 
bond dimensions on are needed, but even for infinite sys- 
tems excellent approximations are possible with moder- 
ate bond dimensions if the ground state fulfills an area 
law for the entanglement entropy |25| . 
For Hamiltonians 



H = 
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a similar ansatz as for the quantum state (pi leads to the 
concept of matrix product operators (MPO) [T8H2T] 
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Many relevant Hamiltonians are represented by MPO 
with relative small bond dimensions. For our purposes, 
it is important to mention that this is often also true for 
Hamiltonians with long-range interaction terms (see, e.g., 
|21|). A recipe for the explicit construction is explained 
in Appendix ^ In the case of a translational invariant 
Hamiltonians (see remark below), the MPO can be built 

are identical 

1. This allows us to drop the index in 
the square brackets except for the leftmost and rightmost 
tensors. 



in such a fashion that all tensors Hfi , lAt 
for 2 < i < n 



II. BASIC ALGORITHM 

In this section we review fundamental concepts |T| and 
the iMPS algorithm as presented in reference [55]. For 
this algorithm to work not only the Hamiltonian but also 
the ground state have to be translational invariant. The 
extension to ground states with broken translational sym- 
metry will be introduced in section |III| 



A. MPS and MPO 

In this paper we deal with spin chains. The quantum 
state of a spin chain is determined by the inner degrees 
of freedom of its components 
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Further more, all tensors H are independent of the total 
number of sites. As a consequence, the MPO of a trans- 
lational invariant Hamiltonian for n sites can easily be 
augmented to n + 1 sites by just inserting another tensor 
H. 

Apart from an efficient representation of quantum 
states and operators, MPS and MPO also provide an 
efficient way to calculate expectation values. For details, 
we refer to the literature (as, e.g., [1]). 

Remark: In this paper we apply the term transla- 
tional invariant also to finite systems (with open bound- 
ary conditions) and their Hamiltonians which are used in 
the iMPS algorithm to approach the infinite case. 
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B. Overview of the iMPS algorithm 



Since the size of the tensor A* 



grows exponentially 



with the number of sites, a more economical representa- 



Any algorithm which deals with infinite MPS could 
be addressed as iMPS algorithm. In this paper we use 
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this term exclusively for algorithms of the type as pre- 
sented in reference [22 . This algorithm aims at finding 
an MPS representation for the ground state of an infi- 
nite one-dimensional quantum system. In its plain ver- 
sion the iMPS algorithm takes translational invariance 
for granted such that all sites behave equally. Thus all 
we need to construct the entire state is a perfect descrip- 
tion of one site and its entanglement features with its 
environment given by the rest of the system. This envi- 
ronment is dominated by nearby neighbor sites while the 
influence of sites far away can be neglected in any one- 
dimensional system with an asymptotic decay of corre- 
lations faster than r _1 . Therefore the environment built 
up by an infinite system can be simulated with a finite 
system. Correspondingly the center site of a sufficiently 
large but finite system provides a good approximation for 
its counterpart in the infinite case. 

The iMPS algorithm is built upon a finite system, 
which is iteratively enlarged by inserting new sites into 
its middle. Since we express quantum states by MPS 
each of these new sites is represented by an individual 
tensor Before a new tensor A[ n ] is inserted it is 

optimized such that the resulting energy 



E = 



is minimized. Hereby all previously inserted tensors 
Au <n i are left unchanged. Of course, these local opti- 
mizations of the new tensors j4[„] are generally not suffi- 
cient to find the lowest energy state of the entire system. 
What we are supposed to get is a ground state approxi- 
mation which might be bad at the outer edges but close 
to the center, it should become better with each new 
tensor inserted. This is all we need to obtain an ade- 
quate description of the center site's environment, since 
the influence of the outer sites fades away with distance, 
anyway. Therefore, we expect the environment of the 
center site to converge towards its infinite counterpart, 
and with that the new tensors inserted into each 
round should converge, too 



A, 



A 



converge 



,!]■ 



(7) 



For ground states which violate translational invari- 
ance, it is no longer given that all sites behave equally. 
At this point, our intuitive argumentation breaks down 
and the algorithm generally fails to converge. We will 
study this case in section For the time being, we 
assume that the algorithm ends up with a converged ten- 



sor A\ 



converged] ■ 



With the help of this one tensor, the 



entire iMPS can be constructed using the rules we will 
encounter in III CI 



1. Long-range interactions 

An important ingredient for a fast computer code is 
a clever book-keeping of the interaction terms, which al- 



lows us to save many calculations due to recycling. In the 
case of long-range interactions, this task becomes tricky, 
since the iMPS algorithm permanently splits the MPS 
and adds new sites. By this means, the distances be- 
tween sites on the left and right halves change each time 
and with them all distance-dependent interaction terms. 
Nonetheless, for translational invariant Hamiltonians re- 
cycling can still be done in a well-arranged fashion by 
encoding the Hamiltonian into an MPO as in equation 
([5]). Every time the MPS is enlarged by a new ten- 
sor, the MPO is enlarged by the standard tensor iJ^/ s Mr , 
too. This simple procedure automatically corrects all 
distance-dependent interaction terms. 



C. Constructing the MPS 

So far, we just mentioned that the iMPS algorithm 
constructs the MPS by constantly inserting new tensors 
in its middle. We will now specify on that. First, it is 
convenient to treat the MPS as divided into two halfs, left 
and right from the center. Each time the optimization 



procedure described in II D provides a new tensor A" lC " r , 
we have to decide into which half A^ l ° r is to be absorbed. 



(6) This is done under the following rules: 



1. According to the half in which A" iar is to be ab- 
sorbed, decompose it (Fig.[TJ(i)) as 
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where the asterisk denotes the complex conjugation 
(see the upcoming equation (10 1 for details). 



2. With each new tensor A overwrite the A of the ten- 
sor before. 

Each A[ n ] is optimized in such a fashion that it compen- 
sates for the overwritten Ar n _ij. 

For the orthogonalization ^ we proceed as follows: 
In a first step, we write the tensor A" l0lr as matrix A m ' a 
where m is a multi-index. If the tensor A" lUr is to be 
absorbed into the left MPS half m — (s, ai) elsewise 
rn = (s,a r ). For the leftmost and rightmost tensors of 
the MPS, which have only two indices, m = s. The in- 
dex a corresponds to the leftover index of A" ,ar , which 
is not in m. Next, the matrix A m ' a is decomposed into 
an orthogonal part Q and a "rest" part A. Different de- 
compositions would fulfill this task, but for the working 
of the algorithm it is best to resort to a singular value 
decomposition 



A = U-D-V^ 

= u -v\-y -d-v\ 

Q A 



(10) 
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Rewritten as tensors we end up with equation (181 . 

A consequent application of these rules yields an MPS 
built of orthogonalized tensors Q and one single matrix 
A from the very last A[ new ] in the center (Fig.[l](M)). 



Q 



l]ai 
[L] Sl 



Q 



[k] antti 



[L]s k 



; • A c 



Q 



[k+l] a k a k 

[R] s k+ i 



! [R]s n 



(11) 

This MPS standard form is numerical robust and has an 
easily calculated norm |)-0|| = \J (ip\ip) ■ To see this let us 
multiply QQ?* ■ ■ ■ QK"*- 1 "' (the left half of equation 



(111) with its complex conjugated 
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and analog for the right half. Thanks to equation Q 
all Q ■ Q*-pairs turn into ^-functions and the norm of 
the MPS (111 equals the remaining ||A|| (see Fig.[I](m)) 



which also equals the norm of the last inserted tensor 
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Figure 1: (i) Diagrammatic representation of the two decom- 
positions of A according to equation (JSJ) . Vertical legs cor- 
respond to physical indices while horizontal legs belong to 
the auxiliary indices. Connected legs are summed over, (ii) 
Resulting structure of the MPS before and after the de- 
composition of the last j4[ new j. (Hi) MPS version of ( , i/>|?/>}, 
where the turned over MPS symbolizes the complex conju- 
gate. Due to the orthogonal decomposition ( 10 1 the tensors 
in the left and right boxes generate the i den tity (12 I, which 
allows an easy control of the MPS norm ( 13 1. 



Effective operator \ 



D. Tensor optimization 

The iMPS algorithm is an iterative procedure. As de 
scribed in section 



II B 



new tensors are constantly 
inserted into the MPS,' which represents the finite state 
ip. Each of these new tensors A^ is optimized such that 

the energy E = ^ctct^ is minimized. As we will discuss 
in more detail below, this can be written as 



(^[n]|U[„]|^4[„] 



(14) 



(^[n]|I[n]|^[n 

\A[ n ]) is the vectorized form of the tensor A[ n ], that 
is |A[ n j) = j4j n j = A?^ with the multi-index i 
(a t ,a r , s). 



Hu is an effective operator built of the MPO rep- 
resentation of the Hamiltonian H (|5| and all MPS 
tensors of (ip\ and except for the two new A\ n y 

I[„] is the identity operation thanks to the orthog- 
onalized standard form of the MPS, see ( 12 1. 



In III CI we mentioned that it is convenient to treat the 
MPS as divided into two halves: left and right from the 
newest tensor Ar n i. For the same reason we decompose 

H[„] into a left half L^ l °" and a right half i?^j MrQr , 
which are connected by a sing le MPO tensor ^ 
corresponding to the new site (see Fig.^l(m)). 



L 
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with i = (ai,a r ,s n ). 

Since the iMPS algorithm is an iterative procedure, 
L[„] and are built up iteratively, as well. Suppose we 
intend to absorb the At n _u of the previous optimization 



step into the left half. First, we use equation 
the orthogonal tensor Qr/i J x ■ With that 
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where the asterisk denotes complex conjugation. In this 
case, where the tensor Ar„_xi is absorbed into the left 



half, the right half stays unchanged i?r„i = R\, t 



Ob- 



verse, if we decide to absorb A^^ into the right half, 



Equation ([hi]) is solved by setting \A[ n ]) equal to the the left half stays unchanged and R becomes 
lowest eigenvector of H[ n ]. 

Remark: We will repeatedly use the notation \T) or R 
(T\ for a vectorized tensor T. 
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Figure 2: (i) Diagrammatic representation of the Hamilto- 
nian MPO (ii) MPS and MPO realization of (ip\H\tp) 

(compare with Fig.[T|. The boxes indicate the left and right 
halfs |l6f and R<»r<* r ^ ^ g ame ob j ect ag 

above with contracted inner indices of L ai ^ iai and R^r^r^r _ 
The object in t he H-shaped box corresponds to the effective 
operator EI 1 15 1. 



E. Algorithm 

After having presented the decisive ingredients of the 
iMPS algorithm, we like to emphasize the steps one actu- 
ally has to perform on the computer. The algorithm con- 
sists of an initializing procedure (see II E 2 1 and an itera- 



tion loop, which is repeated until convergence is reached 



-4, 
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[converged] 



0. 



The tensor A 



converged] 



is all we 



need to calculate expectation values. We do not hold any 
copy of the MPS we are calculating. The only objects 
stored (in the purest version of the algorithm) are the 

actual versions of Lff 1 " 1 (pi, R<^ a - pj| an d A 



aia r 
n] s 



3. Use (|8| and (10 1 to decompose AY, and gain Q\l] 
or Q [R] . 



4. Use (16) or (17] to get L£+i] and R [n+i] ■ 



2. Initialization 

At the beginning we have to initialize the values of 
and R<^ a -. This can be done with the help of 
an exact solution ip for a small system best with an even 
number of sites n = 2k (we mostly used n — 8 or 10 for 
two-level sites). The state ip is split into a left and right 
part, which allows to calculate L and R. In more details, 
note the following: 

1. Write the Hamiltonian of the small system as ma- 
trix with multi-indices _ff"( s i'" s ™)'( Sl "" Sn ) and solve 
for ^(si-'-S"). 

2. Split the multi-index in two multi-indices, giving 
l p(s 1 "-s k ),(sk+i-- •s„) ) wmcn can be interpreted as a 

matrix. 

3. Use a singular value decomposition ip = U ■ D ■ 
(or Takagi's factorization ip = U-D-U T Hip = ip T )- 

4. Interpret the index structure of U as U 1 ^ = JJ l ' ak = 

Jj(si---s k ),a k 



5. L a ^ kOL 

■s' k ,fj, k 
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6. Use V analog to [/ to calculate R a >- a r . 



Comparing with (11 1 we find that 
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III. BROKEN TRANSLATIONAL INVARIANCE 



1 . Loop 

1. Calculate the new ^X"*" = A]^ with i = (ai, a r , s). 
Therefor 



(a) Use (15 1 to calculate 



(b) Set equal to the lowest eigenvector of H^ 1 

2. Decide, whether to absorb A 1 ^ into the left or right 
half (e.g. even steps left, odd steps right) and act 
accordingly in the following two steps. 



In this section, we present some new concepts for the 
iMPS algorithm which arise from the need to deal with 
spontaneously broken translational invariance. A princi- 
pal shortcoming of the basic iMPS algorithm is its failure 
to converge in such cases. To overcome this deficiency, 
we introduce the superposed multi-optimization (SMO) 
method in section |ilI B| Once convergence is restored, we 
turn our attention in section [ill C| to the question of how 
to obtain a specific solution out of the ground state man- 
ifold degenerate due to broken translational invariance. 



In addition, in section III D we treat the special case of a 
non-degenerate ground state which is separated by a very 
small energy gap from a state that breaks translational 
invariance. 
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A. Preliminary considerations 

Just one tensor A [converged] suffices to construct an en- 
tire iMPS. At first sight, one might therefore think that 
such an iMPS is only capable of describing states where 
all sites behave equally, which is no longer true for states 
which break translational invariance. But still, also these 
states can be handled. This is due to the construction 
rules presented in section |fl C These result in an iMPS 
structure given by equation (11), where the matrix A 
marks a special position (and with that breaks transla- 
tional invariance) if it can not be commuted to its neigh- 
bor sites. 

In 



II B 



The real problem is to find Ar convergod ] . 
already mentioned that our argumentation in iavor of 
the convergence At n ] — > A[ con verged] is no longer valid in 
the case of broken translational invariance. The iMPS 
algorithm is grounded on local optimization and therefore 
it is vulnerable to locally altering states, as they appear 
on a physical level for states with broken translational 
invariance. In this case, simple local optimization will 
not result in a global optimal fixed point At, 



converged 



1. Known solutions 

When we write about the breakdown of translational 
invariance, we mean that the state is no longer invariant 
under the shift of one site. Still, the state can maintain 
invariance under the shift of k sites. If k is known, one 
can introduce new super sites where one super site en- 
compasses k of the old sites. Now, the system is transla- 
tional invariant with respect to the shift of one super site. 
This involves that we have to optimize tensors which rep- 
resent the super sites. Due to the exponential increase of 
the physical dimension, this method is practical for very 
small k only. To avoid the scaling problem, Crosswhite 
|22| suggested to use an MPS ansatz for the super sites. 
In practice this means, we insert k old sites at once and 
optimize the corresponding tensors. We are not aware if 
this was ever tested successfully. Aside from that, one 
still needs a prior knowledge of the value k. 

Alternatively, one can extend the standard MPS struc- 
ture with auxiliary tensors [26 , 27J , which allow us to in- 
troduce a symmetry breaking element for the price of a 
non-linear optimization. To our knowledge, this has not 
been tested for long-range interactions, so far. 



B. Superposed multi-optimization (SMO) 

Our solution of the convergence problem induced by lo- 
cally altering states does not depend on any prior knowl- 
edge and stays within the standard MPS framework. The 
key idea is to wash out local dependency of the opti- 
mization by choosing each new tensor Ar n i such that it 
minimizes the sum of the energy of exponentially many 



different MPS instead of just one. These MPS are dif- 
ferent ground state approximations to the qualitatively 
same Hamiltonian applied to systems of different sizes. 
All necessary minimizations can be joined in a superpo- 
sition and solved by one optimization. The time relevant 
steps stay the same as in the single MPS optimizing al- 
gorithm presented so far, so there is no noticeable loss in 
speed. 

As explained in the following, after each optimization 
round, the number of MPS joined in the superposition 
increases by a factor of 4. Thus, in the nth round, the 



optimization ( 14 I gets formally extended to 



f 4 n- 



f 4 n- 



min ^{ipilHilipi) -^min ^ (A [n] |H [n] i\A [n 



4 n- 



= min (A[ n ]\ H[ n ] i\A[ r , 



= min^(A w |H N |^ w )J. (18) 

The MPS representing the \tpi) are of different length and 
the position of the tensor A\ n i is no longer in the center 
but varies from MPS to MPS. In the basic algorithm, the 
tensors Ar n i experience quite individualized environments 
and the optimization adapts to these local circumstances. 
In the modified algorithm, A[ n ] faces exponentially many 
different environments averaging out local effects and em- 
phasizing common global features. This enforces heavily 
the desired convergence Ar n i — > A [converged]- 



1. Modification of the algorithm 

Formally, the superposition Hj„i = X^=i ^[n] i m 
equation (pi) is based on 2 2 ("" 1 ) = 4™" 1 different MPS. 



These MPS are not hand picked, but indirectly gener- 
ated by the algorithm. The only modification needed to 
create such a superposition concerns the left and right 
halves L a '^' a ' and R a 'r^ a - . We still use equation dl6 > 



and ( 17 1 to perform the iteration steps — > L?\ ' 



and R, 



R 



but afterwards we add the new 



and the old results to achieve superpositions 



R 



*~ R \n-1] 



[n] 



(19) 



Since this is a simple addition of two tensors, the struc- 
ture and size of L,^ 1 1 and iir^j r stay the same and 
do not entail any computational complications. 

In the basic algorithm we have to decide at each itera- 
tion step whether to absorb the tensor A[ n _i] into the left 

halfL^'" 1 ((16) or into the right half K^f r<Xr (^7). Only 
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the half of choice is modified. Now, we symmetrize the al- 
gorithm and modify both halves in each step. With this 
modification, the operator sum H[ n j is calculated with 
one single use of equation (15 1. As a further advantage 



The first question can be answered more formally. 



of the symmetrization the modified iMPS algorithm can 
now take advantage of mirror symmetric Hamiltonians. 
As this subject is a bit off-topic, we refer the interested 
reader to Appendix [F] 

Since each of the tensors ^[1], ■ ■ ■ , ^-[n-i] is absorbed 
into the left half L and the right half R the longest MPS 
encoded in the operator superposition of H[„] contains 
2(n— 1) tensors (where we neglect the contribution from 
the initialization routine explained in section |IIE2| and 
the hole for the new tensor ^4[„]). But in the general MPS 
each of the 2(n— 1) tensors only appears with a probabil- 
ity of 50% due to the addition of the new and old L and R 
in equation (19 1. The possible tensor combinations give 
rise to the heralded 2 2 (™ -1 ) different MPS from which 
( 2 (n-i)) are of len g th k - More Precisely Q^) ■ Q^) 
of these MPS contain tensors left and k^ tensors 
right from the hole for the new tensor A< n i. 



2. Comments 

The crucial observation is that the iteration steps to 
perform are always the same independent of the tensor 
position and the size of the MPS. Therefore all the dif- 
ferent MPS can be optimized together combined in a su- 
perposition. The reader who is more familiar with finite 
MPS calculations might wonder about the complete loss 
of information concerning the single MPS in the super- 
position, which comes along with equation ( 19 1. We have 



to remind ourselves that the main objective of the iMPS 
algorithm is to get the tensor A [converged]) which suffices 
to construct the infinite MPS. The finite MPS are just 
tools to obtain this tensor. Once we have it the finite 
MPS are no longer needed. 

Another interesting question is whether the operator 
sum Hr n i in equation (18 I constructed via equation (19) 



is really suitable for a variational ansatz. This question 
has two aspects. 



1. Is every optimal iMPS solution of equation ( 18 ) also 
a minimum of the Hamiltionian? 

2. Does the algorithm always converge to this optimal 
solution? 



The intuitive argumentation given in section |II B| also 
suggests that the answer to the second question should 
be yes, but actually even for the well-established DMRG 
algorithm the answer has to be no, since otherwise NP 
hard problems could be solved. Nonetheless, it is a mat- 
ter of fact that DMRG converges extremely well for most 
practical purposes. The same "practical proof" can be 
given for equation ( 18 ) as demonstrated by our applica- 
tions (section |v| . 



Equation (18 I represents the sum of exponentially many 
energy terms. The lowest conceivable value of this sum 
is reached, if each energy term takes its individual min- 
imal value. If all individual energies are minimized, the 
obvious answer to the first question is yes. Hence, we 
have to ask, whether it is possible to minimize all in- 
dividual energies at once, having in mind that all MPS 
involved are created indirectly via equation (19 1. With 



finite MPS this might only be possible up to a certain rel- 
ative error. But for the limit of infinite MPS this relative 
error shrinks to zero and the problem is trivially solved 
by uniform iMPS, i.e., in the case where all tensors stem 
from the same ^.[converged] ■ Then, all iMPS appearing in 
equation (18) look alike and either none or all of them 
minimize their Hamiltonians. Even in the case of broken 
translational invariance, one can always find at least one 
translational invariant ground state, which can be writ- 
ten as uniform iMPS and hence optimizes all terms in 
the sum at once. 

Next, we like to further inspect the numerical conse- 
quences of equation ( 19 1 for the different MPS which are 
part of the operator sum M.\ n i. The tensor Ar n _;n was 
absorbed into exactly half of the superpositions encoded 
in L[„] and in R\ n -\. Since Lr n i and R^ are the building 

blocks of 
guished: 

1. Ar„_x] was neither absorbed into L[ n ] nor into R[ n ]. 

2. Ar„_xi was only absorbed into Lr n i. 

3. Ar n _ii was only absorbed into R\ n ]. 

4. Ar„_xi was absorbed into both halves L[ n ] and R[ n ]- 

Although we never experienced any practical problems, 
the cases 1 and 4 are, at least from the theoretical point 
of view, a bit troublesome. In case 4 the tensor ^4[„_i] is 
inserted twice. But Ar n _i] was never optimized for dou- 



[„] (15 1, four subsets of H[„] can be distin- 



ble insertion. Close to the end, when At r ,i — > Ar, 



rged] 



is almost achieved, this should pose no problem. Mean- 
while, at an early stage the effect should be more severe. 
On the other hand, even in the basic algorithm, the MPS 
description is not perfect - especially not at the begin- 
ning. 

Case 1 might seem trivial, since everything stays the 
same. Potential difficulties arise in the superposition 
with the other cases. According to equation Q the ten- 
sors A are decomposed into Q and A and only Q is ab- 
sorbed. The matrix A is overwritten with the next A 
(respectively Q). In the basic algorithm, this is easy to 
justify: The next A[ n ] can compensate for Ar n _ii. But, 
a perfect compensation can only be achieved for one A, 
not for many of them. Here is the problem: In case 1, 
old A[„_ 2 ], A[„_3] i ... of the previous steps are conserved, 
while in 2., 3. and 4. a new A[„_i] comes into play. All A 
have to be compensated for. The stronger the A alter, the 
less adequate is their compensation. The variation of the 
A can be reduced by enforcing \\A\ n ] — Ai n _-u | to be small. 
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Towards the end of the optimization ||A[„] — A[ n _i]|| is 
small anyway. At an early stage one might have to resort 
more strongly to the convergence enforcing method we 
will discuss in section IIV A 21 



C. Selecting a specific ground state 

We have seen how to ensure convergence in the case 
of broken translational symmetry. But so far, we have 
no control to which of the degenerate ground states the 
algorithm converges. Some of these ground states might 
be more favorable for our purposes than others, and we 
now answer the question as to how to obtain them. Any 
further degeneration aside from broken translational in- 
variance is excluded from this consideration. 

Let us look at two fully converged MPS A and B 
where B be a representation of the wished-for ground 
state which fits our purposes best, while A stands for 
any ground state to which the algorithm actually has 
converged. According to equation (111, both MPS have 
the following structure 



A _ n°-2 a -i n a -i a o \a ao n»o a i noio 2 
— ■''*[£]«_! "*[i]s ' A ' ^[«]si ' ^[R]s 2 " ' 

n _ a- 2 a-i a-iSo fa a a ai a 1 a 2 Mn-, 

D _ q [L]s^ 1 1[L]s ? Q[R] S1 1[R]s 2 \ AU > 

In Appendix [C] we show that it suffices to replace the 
matrix A a ° Q ° in A by the new matrix j a ° a ° to obtain an 
MPS which represents exactly the same physical state as 
B 



B 



Q 



OL — 2 OL-\ 

L]s-i 



{ [L]so 7 



'[R] si 



In other words, we do not need to take care to which 
ground state the algorithm converges, since after it has 
converged, we are able to transform the obtained solution 
easily into any other. We do not even have to know B, as 
long as we have a description such as, e.g., "the ground 
state with the highest expectation value for the operator 
X". All we have to do is a one-time optimization of the 
new matrix -f aaa ° under the desired side condition. 



1. One tensor update versus multi tensor update 

So far, we focused on uniform MPS which result from 
an algorithm that inserts one new tensor each round. 
Crosswhite [55] suggested that one might also insert a 
certain number of q tensors per round in the form of 
a small MPS. Although we are so far not aware of any 
successful practical applications of this ansatz, it is worth 
having a closer look. For the single-site algorithm to 
work in the presence of broken translational, invariance 
we introduced the SMO method, which washes out local 
variations and thereby fortifies convergence. Still, the 
SMO is a general method and could also be implemented 
in an g-site algorithm. 



As we have seen in the section above, the single-site 
iMPS algorithm will come up with a solution that en- 
codes all possible ground states. This abundance has 
its price. Given the situation that we know the period- 
icity q of the ground state of interest, we could use an 
iMPS algorithm which inserts q sites at once. This would 
enable us to find some specific lowly entangled ground 
states which could be expressed by a non-uniform MPS 
with a far smaller bond dimension. Since q translation- 
ally shifted copies of such a non-uniform MPS always 
allow us to construct a uniform MPS, the maximal gain 
in bond dimension is given by a factor q and the maximal 
difference in the entanglement entropy of the half chain 
is AS = log 2 (g). We could confirm this difference for the 
model studied in the applications (section [V]) varying the 
matrix A (20 1 over the set of ground states, as described 



in the section above. 

From the perspective of the needed bond dimension, 
the multi-tensor update is clearly superior to the-single 
tensor update for systems with a periodicity q > 1. Sill, 
in this paper we favor the single-tensor update, which 
needs no prior knowledge of the periodicity and results 
in a well-converging algorithm, which has proofed its re- 
liability in practical tests. 



2. Degenerate tensor 

In the case of broken translational, invariance one can 
jump from one ground state solution to another just by 
changing the matrix A, which is part of the bigger tensor 
A[ n ] (8l. Hence, different minimize (^4[ n ]|H[„]|A[„]), 
i.e., Ai n -i is degenerate. The iMPS algorithm aims for the 



convergence Ai n i — > A 



converged] • 



Without precautions, 



this convergence might be undermined towards the very 
end by an A[ n ] which jumps from one solution to an- 
other. At a first glance, this does not seem troublesome, 
because all solutions A^ could jump to are good solu- 
tions. Nonetheless, due to imperfect numerics this jump- 
ing might also occur into Ar„i of minor quality. This 
effect is not fatal, but it still might turn an otherwise 
perfect result into a less accurate one. 

To suppress this effect, we can resort to the conver- 
gence enforcing method we will present in IV A 2 In 
addition, we will describe in |IIID| and |IVD | a numeri- 



cal method to enforce a translational invariant solution 
which eliminates the above-mentioned degeneration. 



D. Translational invariant ground states and local 
minima 

In this section, we consider possible convergence prob- 
lems due to translational invariance breaking states 
which lie closely above the non-degenerate ground state 
level. In such cases, the infinite system still provides a 
translational invariant ground state, while for finite sys- 
tems even small alterations of the energy spectrum due to 
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boundary effects suffice to favor a ground state with bro- 
ken translational invariance. Since the iMPS algorithm 
is based on growing finite systems, it might start out 
converging into a false minimum and get trapped there. 
Even if the algorithm escapes out of this trap later, it 
supposably costs many optimization rounds and signifi- 
cantly slows down convergence. 

To avoid these problems we suggest to modify the algo- 
rithm such that it only converges to translational invari- 
ant states. This is no limitation: In the case of an unique 
ground state, the sole solution has to be translational in- 
variant, anyway. If the ground state level is degenerate, 
one of the solutions is translational invariant and accord- 



ing to equation (21 1 we can still transform it into another 



type of solution after the algorithm has converged. 
Whether the fully converged MPS A, 



A _ fl a -2 a -i /l a -i a i) \a a r> a o a i r\ a l a 2 (00\ 
— ' ^[L] s-i ^[L]s ' A ^[R]si '^[R]s 2 " ' \ LA > 

is translational invariant or not depends on its matrix A. 
At this point we should be more precise and write Am 
or \\r\ , depending on whether A stems from a left or a 
right decomposition (]§]). Actually, as a consequence of 
decomposition (|8l the MPS A is translational invariant 
if the left and right version of A are identical 



\l\ — \r\ — A ^ Q[ L ] • A — A\ 



converge- 



d ] — A • Q[r]. (23) 



In this case A can be commuted to any position and hence 
does no longer mark any specific site of the MPS. This is 
what we are aiming for. 

In order to end up with an A [converged] where Art] = 
Arm we alter the minimization routine which computes 
the tensors A such that solutions with small differences 
I|A[l] - X[r]\\ i-e. big overlap (ArnJA^j) are preferred. In 
the long run this should accumulate to A[^j = X[r]- 

As a first straight forward way we tried to extend the 



minimization (18 1 of (A|H|A) to 



min((A\a\A)- 1[X y{\ [L] \\ [R] )) (24) 

with a suitable coupling parameter 7^] . This is no longer 
a simple to solve bilinear problem since one needs to per- 
form the decomposition (10) to get Am and A 



[R]' 



To 



avoid this complication and restore bilinearity we tried 
to resort to the easily calculated approximations Am and 
A[/j] (E2| derived in the appendix |e| but the results we 
obtained in this way were not very convincing. 

In section |IVD| we introduce a less conventional ap- 
proach which turned out to work far more satisfyingly 
for us. Instead of extending the minimization of (A|H|yl) 
by a new term as suggested in equation (24 1, we alter 



the routines of the iterative eigenvector solver we use to 
solve it. The modus operandi of these solvers is reviewed 
in section [IV C| Until after then, we suspend further ex- 
planations. 



IV. ENHANCED ALGORITHM 

The considerations of the last section were mainly con- 
ceptual. The only actual change of the algorithm we per- 
formed is given by equation (19), which incorporates the 
SMO method. In this section, we delve far more into 
numerical details and extend the algorithm by further 
routines to make it more efficient. A reader not inter- 
ested in technical details of the algorithm might proceed 
directly to section IV] 



A. Enforcing convergence 

The goal of the iMPS algorithm is the global conver- 
gence j4r n ] — > ^4 [converged] ■ This property has to emerge 
over the long term, while it is not part of the evaluation 
system of the local minimization from which each A[ n ] 
is drawn. As a consequence, small local improvements 
might be purchased with strong fluctuating Ar n i coun- 
teracting global convergence. In an unstable scenario of 
overcompensation, these fluctuations might even inflate 
in a fatal manner. To prevent this from happening we 
extend the algorithm by two methods. The first method 
(superposition method) aims at attenuating the influence 
of problematic A^ on the ongoing calculations, while the 
second method (gain function method) directly modifies 
the optimization routine such that excessive variation of 
the A< n j are suppressed. Both methods are complemen- 
tary and worked well together in our calculations. 



1. Superposition method 



The first method takes advantage of the fact that the 



Hr„i ( 18 1 of the modified algorithm represent superpo- 
sitions of operators. By decreasing the weight of those 
contributions to the superpositions which contain prob- 
lematic A\ n ] , one can ensure that excessive fluctuation 

of the A[ n -\ do not spread to the level of the H[ n _|_i] and 
with that inhibit a chain of overcompensation. We re- 
mind the reader that At„t is absorbed into L 



and i?[ n+ i] ( 17 1 before equation ( 19 ) is used to build up 
superpositions. This latter equation is now replaced by 



J [n+1] 



(16) 



L [n+1] *~ L [n] + M«] ' L [n+1] 
K [n+1] <~ H [n] + ?[n] ' H [n+1] ' ( 2b > 



The only new ingredient compared to equation ([19]) is 
the adjustable weight 1 ^ £[„] > calculated as 



e W =mm 1,— 



(26) 



where AAu measures the deviation of Ar„i and (AA)[ n j 
is a weighted average of the previous deviations. Each 
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time the deviation AAr n i exceeds the average value 
(A A) r„i , £r n ] gets smaller than 1 and with that the weight 

of all contributions of H which contain Ai n ] is reduced ac- 
cordingly. 

To measure the deviation AA[„] we need to define a ref- 
erence tensor A^j such that Ai[ n ] = ||A[„] — v4|^j fer ' ||. 
In order to avoid unnecessary fluctuation of this reference 
tensor we use the same trick as above and define A[ IC ! fcl ' 



iteratively as a weighted average of the previous ^[o^j 



<n] 



(A^4)[„] (28 1 ensures that £[„] (26) is lower bounded 
around 0.1 . 



Assigning the parameter < q n i ^ 1 ( 30 1 allows us to 



shorten AA! 7 ] to a chosen fraction of the maximal value 

[n] 

AAj 7 j~ '. The price to pay for a C[ n ] < 1 is a lesser energy 
improvement A£^J which is calculated as the difference 
between the energy one gets due to choosing A[ n ] — A^p\ 
instead of just taking A[ n ] = A^ er ^ 



.4 



refer] 
n+1] 



1 

N 



A 



[refer] tf . 
[„] + 4[n] ' A [n] 



(27) 



with N = \\A 



refer] i 



The weights used in equation 



( 27 I are the same as in equations (|25|) and ( 26 I 



For equation (26) to work we still have to define the 
weighted average TjAA) r n i . Various definitions are possi- 
ble. As a heuristic choice, we picked the following one: 



(AA) [n] = i • min(0.9 • (A A) [„_i] + 0.1 • AA [n _ 1] , 



1.02 • (AA)[ n _ 



(28) 



with N = 1 - 0.9™. Obviously the term 1.02 • (AA) [n] 
prevents a too sudden increase of (AA)[ n j by limiting it 
to 2% per round. Without this term we get the clearer 
expression (AA) [n] ~ 2™ =0 0-9 n ~ j -AA^, i.e. older AAy] 
lose each round 10% of their influence in the weighted 
average. 

Finally, we remark that we end up in a deadlock if 
= 0. To prevent this from happening we will intro- 
duce the parameter A max in equation ( 30 ) of the upcom- 
ing subsection. 



2. Gain function method 

The idea of the gain function method is to manipulate 
the minimization procedure of (j4[„j |H[ n ] |A[ ra j) by adding 



a gain function i.e. replacing Hr„i by 



AE 



m _ 



<<i|H[n]l43>-MG 
1 - (1 - C[„] 



cfer] I 



\A 



efcr] \ 



AE 



£■[7=0] 



(31) 



Choosing e.g. a 7 which corresponds to cr ra ] « 0.9 reduces 
AAj 7 | by 10% while the energy improvement AEj 7 | is still 
at 99% of the maximal value AElV^ ■ 



The parameter A max ( 30 1 and the entire superposi- 



tion method are designed to intervene only in case that 
Aj4[„] suddenly increases with ongoing n - otherwise they 
have no effect. The parameter cr n i on the other hand al- 
ways effects the calculation if chosen to be smaller than 
1. Generally the q„] should be chosen in dependence of 

AA^j -0 ' (the bigger AA| 7 j~ ', the smaller C[„] and vice 
versa). Just for orientation (not as exclusive choice) we 
give the value we chose for most of our calculations 



0.1; 0.7 + 0.1-log 10 fA^r 0lN )l . (32) 



With that 0.269 < c [n] sC 0.9 since AvlM 2. This 
formula was found heuristically and worked fine for us, 
although more adequate choices might exist. 

When the iMPS algorithm finally approaches its end 
7 becomes very small and its effect might be overruled 
by numerical imprecision. To prevent this, we recom- 
mend defining a lower limit for 7 above the limit of the 
numerical precision. 



HW=H w -7-|<f rl )(<l ferl | with 7 >0, (29) 



where \A, 



cfer] \ 



is the vectorized version of the reference 



tensor defined in equation (27 1. Let A^ be the result of 
the above optimization. Obviously, bigger values for 7 
favor smaller deviations A^4[ 7 ! = ll^! 7 ! — A[ rl f CI '||. 

In the appendix [B] we show how to approximate 7 ef- 
ficiently such that 



AA 



[7] 



■ AA h=0] A 

Tr»l ' max 



(30) 



where < cr n i ^ 1 and A max are parameters of our 



choice. Limiting AA^A by assigning e.g. A max 



= 10 • 



B. Energy overgrow 

If the average energy per site of an infinite state does 
not equal zero, the total energy of the entire state is 
±00. Of course, we never have to deal with an infi- 
nite value since our numeric is restricted to finite sys- 
tems. Nonetheless, a problem remains. In the long run, 
the numeric value of all the information encoded in the 
tensor H stays more or less the same except for the en- 
ergy, which grows with each new site. The tensor H gets 
more and more ill conditioned since the numeric value 
of the energy overgrows other information and thereby 
reduces the achievable precision. To avoid this prob- 
lem, we advise to subtract each iteration step the energy 
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E[ n ] = (A [n ]|H[ n ]|A[ n ]) from the system, 
ing, we recommend to assign 



Em • I 



Simply speak- 



(33) 



But to be of 
encoded into 



L r J 1 1?' and R? 

[n+l] ['i 



L [n+1] 



fl5). 

tensor H^^ r 
in Appendix 
represents a 
tion we add 



any use, this simple assignment has to be 
•Protr |_ ne building blocks of 

This can be done by modifying the MPO 
used in equations (16 1 and (17 1. As shown 
A| the MPO tensor H£f T has a slot which 
ocal interaction term. To this local interac- 



-E h 



C. Minimization routine and information recycling 

With an increasing number of rounds n the succes- 
sive minimizations of the different (Af n ]|Hr n i|^4r n i) be- 
come more and more similar, which opens the oppor- 
tunity to speed up the minimization recycling informa- 
tion from preceding turns. In order to understand these 
ideas (and also those of section IV D section IV E and 
Appendix IB 1} 



we have to review the principles of the 
iterative eigenvector solvers we use [55] . In the MPS con- 
text, these solvers come with the major advantage that 
Hr n i never has to be constructed explicitly; it suffices to 

be able to assemble Hr n ] \A) for any given \A). Further, we 
do not need to perform the minimization to its very end. 
For the algorithm to work, it suffices to perform a limited 
amount of iterations, such that the resulting \A) might 
not be optimal but still significantly improved. Due to 
information recycling, these improvements accumulate, 
such that the optimal solution emerges in the long run. 

Iterative eigenvector solvers are very well suited for the 
outer eigenvalue spectrum. Already with modest effort 
we can expect to find a good approximation |eo) ~ \Eq) 
for the lowest eigenvector of Hr n j . The central idea is to 
project the problem defined on a huge space of dimension 
N onto a much smaller subspace of dimension k <C N 
and solve it there. For this to work, we have to build up 
iteratively a small set {|2li), ■ ■ ■ , |2tfe)} of k orthonormal 
vectors which enables us to express the minimizing eigen- 
vector |t4.["i ) = \Eq) of Hr n i as a linear combination 



l ) = \E Q ) « |e ) 



k 

£ 

»=i 



|Oi> 



Eq — (Eo\W[ n ]\Eo) 

« a! minlt •<2t i |H w 121,) 

_ [min]f ~[n] [min] 



3 

e . 



(35) 



We need to solve for a ™ , which is obviously the 
minimizing eigenvector for the k x k matrix Sy^ 1 = 
(2li|Hr n i|2tj). A possible measure for the accuracy of the 



approximation ( 34 1 is given by the norm of the residual 
vector 



\r) = (H w -e ) |e ). 



(36) 



As long as ||r|| is too big we have to extend the set 
{|2li), . . . , |2lfe)} iteratively by a further vector |2lfe+i). 
Any form of educated guessing for a suitable new |2lfc+i) 
is allowed. The Lanczos [25] and Arnoldi [30 algorithm 
use a different way of calculation but end up with 



|2l 



fc+i/ 



(37) 



with |2lfc+i)_L|2li^j^fc) by construction. 

In contrast to the basic iMPS algorithm, which sets 
|2ti) equal to the lowest eigenvector |er n _i]o) = |^[n-i]) 
of the last round 122 . we choose 



= \A 



cferl 



(38) 



with A[ r< j fcl ' defined in (27 1. This small change allows 



an easy implementation of the method presented in Ap- 



pendix B 1 and should help to improve global conver- 
gence. Both versions are straight forward examples of 
information recycling since |j4[„_i]) as well as |A|^ fcl ') 
are already good approximations for \Eq) (=\E^ )). In 
many cases, we could obtain a considerable speed up ex- 
tending this idea to a few more than just the first vector 
of the set 



|2l 



1+3 1 



A [roforl \ i 



(39) 



1121 



1+3 \ 



refer] \ 

n-j] I 



Further, we observe that all \A^ cic ^) d27k are derived 
from the best eigenvectors of the previous rounds. As 
an additional extension we also tried to include the next 
best eigenvectors \e\ n -i}j>o) of the last round 



|2l 



1+3/ 



(40) 



The improvements we achieved in this way were relatively 
poor. A much more promising way to take advantage of 
the |er n _m) is to use them for an efficient approximation 

/ ~ \ -l 

(JD6|), which 



(34) of the inverse operator D = 



eo 



allows a handy implementation resembling the Davidson 
(or Jacobi-Davidson) [21] method. As a result the update 
equation ( |37[ ) is replaced by the more appropriate ansatz 
( D7 1 . Details are explained in the Appendix [D] 

At the end of this section, we like to caution the reader 
that the methods presented here might counteract the 
methods presented in section |IV A| Global convergence 
and improved local minimization often go hand in hand, 
but not always. If the algorithms indicates to run un- 
stable, one should consider to partially switch off the im- 
provements just presented. This is likely to happen if | eo) 
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is degenerate. In this case, the recycled knowledge from 
the past strongly increases the probability that already a 
shallow optimization suffices to find alternative solutions, 
which might result in unwanted fluctuations as, e.g., de- 
scribed in |III C 2") Usually, this problem is announced in 
advance. In the Davidson implementation, one should 
not resort to eigenvectors with eigenvalues too close to 
the best. Similarly, once the small set of recycled initial 
values |A ^ 6 ^) suffices to get a second best eigenvalue 
very close to the best, it might be wise to abandon this 
method and only use jAf 1 ^ ^) alone. 



D. Enforcing translational invariant ground states 

In this section, we demonstrate the algorithmic real- 
ization of the considerations put forward in section |IIID| 
There, we argued that it is beneficial to push the al- 
gorithm towards translational invariant iMPS solutions 
to avoid getting trapped in local minima. We further 
showed that translational invariance is assured if the 
decomposition 
A 



of the tensor A f conV orgod] results in 



= A[^j] (23 1. This is what we are aiming for. 
The approach we are about to present is not very in- 
tuitive. Therefore, we start our explanations with an in- 
termediate step and introduce a less practical but easier 
to understand procedure which consists of the following 
steps and has to be performed with each new tensor A 
after it has been optimized: 

1. Decompose A into Q[ L ] ■ A[ L ] = A = A[ fl j • Q[ R ] 

2. Define A [sym] = \ (A [L] + X [R] ). 

3. Set A 4r- g (Q[L] ■ \sym] + \sym] ' Q[R})- 

4. Goto 1. 

Due to line 2. this procedure converges towards a tensor 
A with A [L] = A [R] . Further we expect ^[ initial l « A^ nal ^ 



if already A|™ tla1 ' ~ ^[r]"'" 11 ■ Nonetheless, the changes in 
A might be too pronounced to be acceptable. To soften 
this approach one can ignore line 4. and just go through 
1. to 3. once. After that we generally still have A[^j ^ X[r\ 
but with a reduced distance |j Arn — A[/j] || compared to the 
initial value. This is all we need to achieve A^j = A[^j in 
the long run. But the new A is still likely not to qualify 
for the optimizing tensor we are looking for. 

Now, we come to the procedure we really use. Instead 
of symmetrizing the tensor A after its optimization we 
integrate the symmetrization into the optimization rou- 
tine. As recapitulated in section [IV C| the optimization 
routine expresses the vectorized tensor Ar n i = \A) as a 
linear combination 



initial] 



\a) = |st*> • ai 



(41) 



of a small set of basis vectors |2l,;) (34 1. The idea is to 
alter these basis vectors |2li) such that we have a similar 



effect as the procedure above. At the stage of the op- 
timization the Q[^/_r] are s ^ m unknown and we have to 
approximate them by their precursors Qr^/OT- 

The |2li) are created iteratively. In each iteration step 
we first create a new |2li) as we used to do (37 1, ( |D8| and 
then alter it. Therefor we introduce \QL) defined as 



|Sti) 



1 



4 ( Q W 



[n-l] 



A[«]i +X[ L ]i -Q [R] 



n-l] 



with a^h^'OK 1 * 73 

[n-l]*a-y igi \ 7 /3 



Af?l^ = g 



(42) 



where we tensorized the vector |2li) in lines 2 and 3. With 
that, we replace |2li) by an orthonormal version of |2li): 



121, 



121,) <- 



2L; 



■121,). 



(43) 



For a better understanding we insert the |2lj) in the linear 
combination (41 1. As shown in Appendix pi we get 



A 



[sym] 



1 21 
1 

2 
1 

2 



(Qm 

(A[L] - 



+ A 



[syn 



{[R]j 



with 



\r]) 



(44) 



which mimics the effect of the procedure presented above. 
But in contrast to the procedure above the story does 
not end here. The important point to notice is that the 
algorithm can still adopt to the alteration ( 43 ) of the ba- 
sis vectors |2U) and come up with alternative solutions. 
More favorable weights a; than those used in equation 
(44 1 are presumably to be found. Even the |2li) them- 



selves are likely to be different since they are calculated 
iteratively according to the needs of the minimization. 
While there are still enough resources to compensate suf- 
ficiently for the negative effects of the enforced alteration, 
the positive effects should survive since the arguments in 
their favor are largely independent of the a^, |2ti) cho- 
sen by the optimization routine. Still, this alteration is 
a trade off, but we have good reasons to believe that we 
gain more than we sacrifice. 

For practical applications, we only need a few lines of 
code to implement equation (42 1, which is also easy to 
turn off for systems where it is not needed, i.e., when the 
unaltered algorithm shows no tendency to run the risk of 
being trapped in a local minimum. In such a case, the al- 
teration is likely to slow down the algorithm slightly. For 
the applications tested by us the loss in performance was 
only marginal. On the other hand we also encountered 
many cases where the altered algorithm clearly outper- 
formed the unaltered one, which was partially even un- 
able to find the correct ground state within the observed 
run time. 
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Although we strongly recommend to implement the al- solved (see IV C I . To formalize this method, let us intro- 



teration ( 43 1 as presented, one could also use a compro- 
mise and only alter the first basis vector |2li) = A^ fer ', 
which has already a strong impact on the outcome of 
the optimization. This reduced version does not come 
with the need to program a new eigenvector solver. Each 
solver which accepts an initial vector |2ti) will do. In any 
case, the gain function in equation (29 1 is understood to 
change accordingly to the alteration of |2li) = Afe fer '. 



E. Length of the MPS 

After n optimization steps, even the longest MPS en- 
coded in the superposition created by the SMO method 
does not surpass the length I = 2 • n + Iq (where Iq is 
the initial length (see IIE2|). For some systems with 



long-range correlations this might be too short unless n 
reaches some considerably high number, which would go 
along with an extended calculation time. To shorten this 
calculation time two methods might be of help: 



1. Use a tensor A c , 



with small bond dimension 



X [small] until a certain MPS length is reached, then 
increase the bond dimension to its final value X[big]- 

2. Use fast Krylov subspace methods [25] to insert the 
same tensor many times (e.g. 10 5 ) into the MPS. 

A simple and comfortable way to increase the bond di- 
mension from X[smaii] to X[big] is to use an isometric 



X [small] X X[big] 



-matrix u a P with 



■ {u T f a ' = 5 aa ' 



(45) 



and proceed as follows after Ar n i has been optimized but 
still not been inserted into Lr n i (16) and Ri n ] (17): 



[n] [n] 



A^ 4- A?'° r 



(46) 



Next, the new tensor A? l ? r is inserted into L { l ? 1 and 

[n\ s [7i\ 



as usual but without the superposition building 



R 

step (llOj) respectively (|25| . To avoid trapping into a local 



minimum one might also consider to add a small amount 
of noise to A^"/ before applying equation ( 46 ) . 

A possible strategy for the small bond dimension 
X[smaii] is to proceed until convergence has been reached 



A,„i -> A 



converged] ! 



but before the bond dimension is 



increased, many more copies of A> 



converged] 



are inserted 



into the MPS without any further optimization. These 
insertions can be done in the standard fashion or gener- 
ally much faster by projecting the problem onto a small 
subspace, similar to the way the eigenvector problem is 



duce the operator X which inserts one copy of A [converged] 
into Lr„i ((161), i.e., 



1- ■ L\~i — L\ 



n+l}- 



(47) 



With that we build up the Krylov subspace IC r , 
JC r = span {L[„] , 1 ■ L[„] , I 2 ■ L[ n ] , . . . , I r_1 • i[„] } (48) 



and similar with R\ n \ (17). As in 



IV C 



we create an or- 



thonormalized system of basis vectors [£&) 
/ fe-i \ 



1 



i=0 



and calculate 3^ the subspace projection of X 



J, 



(sum 



(49) 



(50) 



Keeping in mind that the subspace projection of Lr n i is 
simply given by the vector lj = ( 1 . . . ) T we find 



L[n+p] = Z P ■ L[n] ~ (p P )ij ■ IjJ ■ 



10 
■h) £, 



(51) 



The number of basis vectors |£j) should be chosen such 
that this approximation is perfect within computer pre- 
cision. Further errors are introduced by an imperfectly 
converged A?' ar ,, and from the energy overgrow ef- 

to conv erged | s toJ to 

feet described in |IV B| which should rule out attempts to 
go for p — ► oo. Still, a small amount of the last two er- 
rors is acceptable since they have a similar effect as the 
aforc-mcntioncd extra noise to avoid local minima. 



V. APPLICATIONS 

In the last section we presented various methods to 
improve the performance of the iMPS algorithm with 
long-range interactions. The main subject was to ensure 
convergence, where special attention was paid to broken 
translational invariance. This so far troublesome case 
can now be tackled mainly due to the newly introduced 
method of superposed multi-optimization (SMO). 

The modified iMPS algorithm has superior conver- 
gence properties compared to the basic version but it 
does not surpass its precision, which is determined by 
MPS and MPO inherited limitations. In cases where 
both versions converge, the quality of the results is iden- 
tical. Readers who are interested in the achievable preci- 
sion of the iMPS method in comparison with analytical 
solutions are therefore referred to the literature [2"2l |2"3"] . 

We checked our algorithm with different models. For 
a reliable basic benchmark, we examined, e.g., states 
with long-range chiral order in the next nearest neigh- 
bor Heisenberg model and found the expected agreement 
with the results given in the references |31[ 132] . 

Here, we will present results for a model of polar bosons 
described by a Bose-Hubbard like Hamiltonian with a 
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Figure 3: The densities p of polar bosons (left) and the corresponding phase diagrams (right) for the ground states of the 
Hamiltonian (531 for U — > oo plotted over t and fi in units of V — 1. Each plot consists of 66049 data points calculated with 



the bond dimensions xmps = 32 for (i) & (ii) and Xmps = 64 for (Hi). The fractions associated to selected phases of the phase 
diagrams denote their p/q values (see main text). 



long-range interaction term. In the thermodynamic limit 
the ground state of this model exhibits symmetry break- 
ing crystalline phases as well as incommensurate phases 
with algebraically decaying long-range correlations. 

The long-range interaction of the Hamiltonian we con- 
sider decays as ~ r~ 3 . To model this interaction with an 
MPO the decay is approximated as weighted sum of 20 



exponential functions 



(r)- 3 *^ ****" 1 r = l,2 > 3, 



(52) 



! = 1 



(see also equation (A17l in the appendix and refer- 
ence [TT?]). 
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Figure 4: The densities p of polar bosons (left) and the corresponding phase diagram (right) for the ground states of the 
Hamiltonian (53 1 with U = V = 1 plotted over t and /j,. Each plot consists of 62194 data points calculated with the bond 



dimensions xmps = 32. The bracketed numbers associated to selected phases of the phase diagram describe their periodically 
repeated occupation pattern (see main text). 
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Figure 5: The blue circles (color online) mark the values of the correlation function (hon r ) for U = V = 1, n = 1.55 and 
t = 0.044 extracted from an MPS with bond dimension \ — 256. For odd distances r — 2j + 1 (right picture) (hohr) is 
suppressed by a factor of roughly 100 compared to even distances r = 2j (left picture). To demonstrate the algebraic decay of 
(hon r ) the function f(r) = a ■ r -0 ' 5 + (non r ^oo) with an adequate a is included on top of both plots as a red line. 



A. Bose-Hubbard model with long-range 
interaction 

We study the thermodynamic limit ground states of 
polar bosons in an one-dimensional optical lattice de- 
scribed by the following effective Hamiltonian |111 155] 

H = V ■ J2 f. \ ]3 ■ ftfc • % + v • J2 h * • ("J - *) 
-H-Y^nj-t-Y^ (c] ■ c j+1 + cj ■ c] +1 ) , (53) 

3 3 

where cj and Cj are the creation and annihilation op- 
erators for a boson on site j and fij = cj • tj. This 



model is characterized by a hopping amplitude t, an on- 
site interaction energy U, a chemical potential /i, and 
a long-range dipole-dipole coupling V/r 3 . For \i > 0, 
the chemical potential favors as many bosons as pos- 
sible in the ground state, while the dipole-dipole cou- 
pling together with the on-site interaction try to avoid 
two bosons coming too close to each other. For certain 
parameter regimes, this interplay allows for translational 
invariance breaking crystalline phases with optimized dis- 
tances between the bosons where q sites accommodate 
exactly p bosons. We will refer to them as p/q-phases. 
The model is known to host an entire Devil's Staircase 
of crystalline phases for t = if the joined potential of 
on-site interaction and dipole-dipole coupling is convex 

urns]. 



1G 



login <Co £ r> 




Log-log plot of the correlation function (c c r ) for 
1, fj, = 1.55 and t = 0.044 extracted from MPS with 



Figure 6: 
U = V = 

different bond dimensions \- Due to the spatial order 
splits into two branches, each decaying oc r~ a + const x , with 
const x < for the lower branch. This constant is a purely 
numerical effect, as can be seen by the scaling with the bond 
dimension \. 



In the following, we will investigate two qualitative dif- 
ferent regimes of this model: U — > oo and U = V. 



1. Devils' s Staircase for U — > cm 

For U — > oo, each site can accommodate at most 
one boson and the effective dimension d e s of the local 
Hilbertspaces reduces to d e g = 2. Figure [3] displays the 
average ground state densities of the bosons and the lo- 
calization of the corresponding p/g-phases. To determine 
the periodicity q of the phases we counted the number of 
eigenvectors of the transfer matrix Tjj?, 



T, 



Q 



[L) t 



(54) 



with an absolute eigenvalue of one. Once q is known p 
follows from the average density. Figure [3] (Hi) shows a 
magnification of the area between the 1 /4-phase and the 
1/5-phase. The biggest phase between these two phases 
is the 2/9-phase, which can be understood as primary 
compromise (2/9 = [1 + 1] / [4 + 5]). In the same fashion 
we find e.g. that the biggest phase between the phase 2/9 
and 1/5 is the 3/14-phase. The maximal detectable value 
of q is given by the bond dimension \ of the MPS, which 
is 64 in case of Figure [3] (Hi). However, since the range of 
\x covered by the different phases diminishes with grow- 
ing value of q most phases beyond q = 30 escaped our 
resolution. The highest value we hit was p/q — 11/52. 
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Figure 7: The Luttinger liquid parameter K in dependence 
of the hopping amplitude t for U = V = 1 and /i = 1.55. The 
parameter K was obtained in two different ways by fitting the 
first order term approximation (551 to the numerical values of 
(hoh r ) as well as (c c r ), both calculated with an MPS bond 
dimension \ = 256. The fits were performed in the interval 
r = 10 ... 60 with only two free parameters, namely const^ 
and K of equation (1551) . 



2. Devils 's Staircase for U = V - 



For sufficient small U the ground states of the Hamil- 



tionian ( 53 1 might accommodate more than one boson 
per site, which allows for new types of Devil's Staircases. 
An example is given by Figure [4] which shows the den- 
sities and phases for U = V = 1. Here, simple trans- 
lational invariance is broken by an underlying occupa- 
tion pattern given by . . . , Xi, 0, 0, 0, . . . with 
Xj = 1 or 2. Of course, for any non-zero hopping am- 
plitude t > we expect fluctuations around this pattern 
such that a more accurate description might be given 
by . . . ,Xi, £j+i) x i+2i £i+3i • • • i which we need in the next 
subsection |V A 3[ At a certain point, these fluctuations 
will become so strong that the underlying pattern is de- 
stroyed, but this is not the case for the entire region of 
Figure [4] 

In the lobes of the new Devil's Staircase the sublat- 
tice . . . , Xi, Xi+2, £1+4, . . . crystallizes in regular pattern 
of single and double occupied sites. These lobes exhibit 
an approximate symmetry under the exchange of single 
and double occupied sites. This is a nontrivial symme- 
try in contrast to the exact particle hole symmetry of 
Figure [3] (z). 

Outside the crystalline phases Burnell [TT] predicted a 
supersolid like phase. In the following we present numer- 
ical evidence which supports this claim. 



3. Supersolids for U = V = 1 

A supersolid is characterized as spatially ordered 
phase which also exhibits superfluid properties. We al- 
ready mentioned the spatial order, belonging to Fig- 
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Figure 8: The entanglement entropy of the half chain S x ex- 
tracted from MPS with different bond dimension \ is plotted 
in dependence of the hopping amplitude t for U = V = 1 and 
fj, = 1.55. The different curves are roughly equidistant, as 
predicted by equation (|56|). 



ure [4] which is given by the occupation pattern 
. . . , Xi, £j+i, £i+2, £i+3i ■ ■ ■ ■ In our numerical studies we 
consider translational invariant superpositions of the 
ground states, where the occupation pattern is still visi- 
ble in the two-point correlation functions as (fiQn r ) and 
(cjc r ) shown in Figure [i] and Figure [s] Due to the spa- 
tial order both correlation functions are split into two 
branches, where correlations belonging to odd distances 
r = 2j + 1 are strongly suppressed compared to cor- 
relations belonging to even distances r = 2j. Further- 
more, both branches exhibit an algebraic decay of the 
same power. For a superfluide, the power of the decay 
of (hon r ) is supposed to be reciprocal to the power of 
(cjc r ). Luttinger liquid theory predicts |35| 



(hon r ) = {fiQhr^oo) + consti • cos (2irpor) ■ r 
(CqC,,) = const 2 ■ r 2K + . . . 



-2K 



(55) 



with po = (hi) ■ These relations are confirmed by our 
numerical results, as displayed in Figure [7j Both cor- 
relation functions give rise to the same Luttinger liquid 
parameter K within a small error range. 

In the thermodynamical limit algebraically decaying 
correlation functions go hand in hand with an infinite 
entanglement entropy of the half chain S, which can not 
be represented by any MPS with finite bond dimension \- 
Nonetheless, it was shown [361 13"7] that the numerically 
obtained entropy S x for such critical phases shows a pre- 
dictable scaling as function of the MPS bond dimension 



$X2 $Xl 




log 2 



(56) 



where c represents the central charge. A demonstration 
of the scaling behavior is given in Figure [8] From this 
numerical sample one obtains ^=256 — S x= iq = (0.218 ± 



drawn from equation ( 56 1 for c = 1 



VI. CONCLUSION 

We have presented several extensions to the basic iMPS 
algorithm for systems with long-range interactions, some 
of them with the potential to be useful in a much broader 
context. A special focus was set on problems arising 
from broken translational invariance. Here, convergence 
was ensured by various means, but mainly due to the 
SMO method which irons out local variations by opti- 
mizing exponentially many MPS simultaneously. The 
new algorithm was successfully applied to calculate de- 
tailed Devil's Staircases and phase diagrams for polar 



bosons ( 53 1 and was also suitable to verify supersolid 
properties. Theoretical restraints as considered in the 
comments |IIIB 2| seem to have negligible influence on 
practical applications such that the new version of the 
iMPS algorithm is a genuine improvement in the sense 
that it can all the old version could plus more. 
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Appendix A: MPO representation for Hamiltonians 

For self-consistency we give a short account based on 
some examples how to construct an MPO representation 
for a given Hamiltonian (see also |18!l2T] ). 

For finite systems with open boundaries the Hamilto- 
nian can be written Q 



SlS2-"«n s[s 1 s' 2 s 2 S3S3 

si ^ S-n—1 s' s n 



(Al) 



First, we need a neat way to write down the explicit 
form of the fourth order tensors H^'^ r . We write them 
as matrices whose entries are matrices, too 



Try-IP- 
s' s 



(H s s ) . (A2) 
As an example we consider the Ising Hamiltonian 



(A3) 
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As we will see below, a possible choice for all in 



equation ( Al I with 2 < k < n — 1 is 



I s s 


s s 


s s 




QS'S 


QS'S 




4 s 





and 1 are vectors over matrices 



l]Ml = / 



-(To; 1 , 



(A4) 



(A5) 



In order to get a better understanding we look at the 
tensor product of the first k tensors. Below (A7) we 
show by induction that 



H [i...k] 



JjlM Ml 
s', si 



tt[2] MlM2 



Tr[k] Mfc-lMfc 



EK-L 
i=l ° 



k-lji) Ji+l] 

V 2 



(A6) 



The resulting vector H I 1 -"*] can be seen as an object with 
three "slots" in which all the relevant information about 
the first k sites are stored. Of course, the number of slots 
corresponds to the bond dimension of the MPO. The first 
slot contains all interaction terms between the first k sites 
only and local terms. Since the fcth site also interacts 
with the k + 1th site, the second slot of the vector passes 

\k] 

on <j l z and finally, the third slot preserves the identity 
I = P'i Sl (g) I S 2 S2 (8 •.. igi F'fc Sfe . For ifW this description 

is easily checked. The tensor W> = IH SS ) (|A4| is 

designed such that it performs the correct induction step 



H [i...k] = H [\...k-i] . H [k] 

-( 



l s k s k Q s k s k QS k S k 

~ai kSk s '" Sk s '* Sk 
-<Jx kSk ai kSk I 3 '* 3 * 



H 



The final tensor -ffj?, Ms given by 



(A7) 



We described the vector fll 1 — fe ] (A6l as an object 
which contains all relevant information of the sites 1 . . . k. 
This description is true not only for Ising interaction. For 
any Hamiltonian we have to identify what these relevant 
information are and design the vector accordingly. As 
convention, we use the first slot of the vector to store 
H [k] the sum of all interaction terms between the first k 
sites only, including local terms. In the last slot we pass 
on the identity. The slots in between are needed for in- 
teraction terms which involve (at least) one of the first k 
sites and (at least) one of the other sites k + 1 . . . n. In 
the case of an Heisenberg chain 



n-l 



n-l 



n-l 



(A10) 



i=l i=l 

the vector H^'" k ' needs five slots 



H [\...k] 



(#1 



[k] [k] [k] 



(All) 



Further slots might be necessary if we do not restrict 
ourselves to nearest neighbor interactions. 

Once we have identified the structure of H^— ,k \ it is 
straight forward to write the first tensor ifW in vector 
form. The matrix structure of all the following tensors 
is constructed column- wise such that the induction 



#-[1-3] = #[l-3-l] . H U] 



(A12) 



is accomplished as in equation ( A7 ) or equation ( A9 1 for 
the final tensor H^ n \ According to our convention, local 
terms, as needed in |IV B[ are always represented in the 
bottom left entry of the matrices. 

For long-range interaction, the recipe given so far be- 
comes problematic. The longer the range of the interac- 
tion, the more information has to be stored in the vector 
#[i— which usually requires more and more slots. But 
there are some exceptions (see e.g. [21 ). An exponen- 
tially decaying interaction needs only one slot - even for 
infinite range. As an example, we look at the toy Hamil- 
tonian 



h: 



<7z 



Sr 1 = ( 

With that we get 

#[l—n] _ jj[l...n-l] , jj 



s n s n 

CTx 



s' s„ 



H ln-1] _ a [n-l] . a [n] _ M 



(A8) 



i=l 3=1 

where i — j — 1 is the exponent of A and not an index. 
First, we have to identify the structure of H^'" k ' 

ff M1 = (wW, Ej^A*-'-^, i). 



= E^ +1I -E^ ] - 



(A9) The crucial observation is, that Ej=i ^ & J ' cr ^ i can 

generated iteratively. The following tensors fulfill this 
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task 



( O s 'i Sl , ai lSl , F'i Sl ) 

J-ai kSh A-F* Sfc O s '« Sk 
o{ kSk l s '« Sk 



\ Q s k s k 
( I s " 



J • (Tz 



s ^ 



(A14) 
(A15) 
(A16) 



In order to encode the polynomial decay as (r)" 3 into 
an MPO we resorted to an approximation as a weighted 
sum of N cxp different exponential terms, i.e. 



(r)" 



(A17) 



In the appendix of [TH] it is shown how to calculate the 
optimal a,i and Aj . For the quality of the approximation 
in dependence of A^xp see [22J . 



Appendix B: Gain function 

We like to estimate the influence of 7 in 

fifr] = H - 7 • |A[ rcfor l) (A[ rcf01 'l I with 7 > 0. (Bl) 

on AW which is supposed to minimize (^Mlff 7 '^ 7 ]). 
For our numerical purpose the following simple approxi- 
mation suffices 



with p[ rofcr l|| = ||S|| = 1 and A[ rcfor ] _L B = const. The 
vector B is extracted from A^ =0 ^, which has to be cal- 
culated first. This might seem inefficient, since we have 
to minimize (^4 1 HT^] \A) twice - once for 7 = and once 
for the final value of 7. The solution is to use a mini- 
mization routine which projects the minimization onto a 
small subspace, as explained in section |IV C| This pro- 
jection has to be done only once but can be used twice. 
Since the projection is the most time consuming part, the 
double calculation is done quite cheap. More to this at 
the end of this subsection. 

Now, we calculate the pseudo energy E using the equa- 
tions pll and fB2]}. 



E = (AW|hW|AW) 

= (l - e 2 ) • ((A [roforl |H|A [rofcrl ) - 7) + e 2 ■ (B\W\B) 
+ 2 • e • y/l-e 2 Re(Al rcfer l \M\B). (B3) 



In the following we approximate e • Vl — £ 2 ~ £■ This 
approximation is not needed but it keeps the calculations 
clear. In addition, the formula we will derive from the 
approximated version is numerically more stable. In our 



program we used the exact version (which we will not 
derive here) only if £(7 = 0) > 0.01. 

With £ ■ y/T^l 3 we, E is a parabola in e. Assuming 
that the apex is the minimum one gets 

. - Re(A^\M\B) 

£ min ( 7 ) = ^ ± . (B4) 

(B\M\B) - (ANfcr]|iHpNfcr]) + 7 



£min for £ m i n <C 1, WC choOSC 



Since ||A[ 7 ! -Afi ei] \\ 

^ II [n] [n] II 

£ m i n in accordance with equation ( 30 1 to be 



£min(7) = min (C ' £min(7 = 0), A max ) . (B5) 

From that 7 is calculated to be 

7 = max (7 [c ],7[ Amax] ) with 
1 - c 



((B\M\B) - (A^ cic ^\U\A^ 



7[A max] = (A^\U\A^) - (B\W\B) 

— • Re(A[ rofer ]|H|B). (B6) 



A^ = v/i-£ 2 (7) ■ A^ + £(7) • B, (B2) m 



A, 



1. Subspace projection and 7 

As mentioned, we have to solve min(A[ 7 l |HW J^H) 
twice: first for 7 = and after that for the final value of 7. 
The idea is to reuse the information gathered in the first 
minimization for the second run. As described in |IV C[ 
the minimization problem is projected onto a subspace. 
The first basis vector of this subspace is |2li) = |A[ refer l) 
Hence the only element of the subspace matrix Sj 
which has to be adopted is f)^! <— fj^i — 7. 
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So, the update of the subspace matrix $j is extremely 
simple to perform but one might wonder whether the 
associated basis vectors |2lj) arc still optimal. For the 
Lanczos [25] and Arnoldi [30 algorithms, which are built 
on pure Krylov spaces, it turns out that the influence of 
7 on the basis vectors gets extinguished, while for the 
extended algorithm presented in |IV C| the optimal choice 
of the basis vectors shows a slight dependence on 7. Nu- 
merically, this is not a serious problem. Nonetheless, 
since the second optimization is the important one we 
recommend to use an approximated value of the final 7 
to construct the basis vectors in the first run. A simple 
and effective way is to use 7[ n -i] from the last tensor 
optimization assuming 7[„i s» 7f n _i]. Alternatively, one 
can solve the intermediate subspace matrices f) and use 
these intermediate results to calculate approximations of 
7 as described above. 



Appendix C: Transformation proof for degenerate 
ground states 

Let A and B represent two different ground states to 
the same Hamiltonian with a g times degenerate ground 
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state level due to broken translational invariance 

A _ ^Q_ 2 a-i «a-iao \a a /)0( 

1/4 — "^[ijs-i "*[i]«o V [-B 
t? _ a-2a-l a_i5 f 5oQ a ai aia 2 /r-<i \ 



naiQ2 
V [-R] s 2 



If degenerations due to further symmetries are involved, 
A and B are supposed to have the same characteristic 
values for these symmetries. This allows us to operate as 
if no further symmetries exist, since all operations we are 
about to use leave these characteristic values unchanged. 
Under this condition we will prove the existence of a ma- 
trix 'y a ° a o such that the ground state B can be expressed 
using the tensors Q stemming from the iMPS A 



Although not an identity, the operator (fj*..^ • (l*f.* Sfe 
acts like such if it is applied from the left on any MPS 
with the structure ^" fc , • , , where 

represents an arbitrary right side of the MPS. Using the 



identity ( C6 1 we get 



i at 



e*fc 

Sfc+l. 



In the following we assume that the two MPS 



(C7) 



B 



Q_2CK-1 

L]s_i '^[L]s 



< ! <*o & ota , 



'[fi] si 



Qai«2 
[R] S2 



(C2) 



where 7 5 ° Q o = u^P . ^<5 . v Sa _ 

Let us start by surveying the elements of the proof. In 
order to show the claimed equation ( C2 1 we will actually 
prove the gauge transformation 



a_2a_i ct — iao 
' ' ' 1[L] s_i • Q[L] s 

Sl • 1[R] S2 ' 



Q 



Q-2Q-1 

[L]s-i 



Q 



[R]si 



L] s 

OL\OL<2 
[R]S 2 



(C3) 



This gauge transformation will be proven for the case 
that the two underlying MPS represent the same phys- 
ical state - which is not given for A and B (Cll. In 



order to use the gauge proof for our purpose we need to 
find one physical state described by two different MPS 
where the first MPS be constructed using the Q of A 
and the second by the q of B. This one state which al- 
lows us to complete our proof is the translational invari- 
ant ground state. According to our preliminary remarks 
this state is unique. Hence, if we succeed to construct 
two different MPS which represent a translational invari- 
ant ground state, we know that they represent the same 
physical state as demanded. We will prove the following 
construction for this state 



a k 



^ Sfc + l... S„ 



a k 

Sl'-'Sfc 



lvJl Sfc + l...S„ 



(C8) 



represent the same physical state. Now, let us apply 

Ifl—Sk 



the operator Qg*... Sk ■ Qti»-s k on this equation. Due 
to equation (C7) the left side stays unchanged. Hence 
the physical state represented by the right side does not 
change either 



) a k 
si—s k 



,l3a k _ 



\*> _ . 
■St— 8k 



re** 

• (J 3 * "- ■ fc* , , 

Sfc + l ■■■S n I 



(C9) 



with ui 
verse 

sured: if the rank of 



1 with M" fc 

s k+i- 



■Mr 1 * 

Sfc + l 



l" k 



The existence of the in- 
= 5 akl3k can be as- 
is smaller than the bond 



dimension ak, the value of ak can be reduced, since in 
that case it turns out to be unnecessarily high. Applying 
Sfl -1 from the right on equation (C9l we end up with 



a k 

si— s* 



Pk 

Sl—Sk 



to 



P k a k 



(CIO) 



T r\Ot,-2<X-\ 



'[L]s_ 



! [L]s ' T 



![R] Sl ' ^[R]s 2 



a- 2 a-i a_itto na a a ai aia 2 
'«[!.] «-i ' »[i]ao ' ' 9 [«]si ■ 9 [fi]s 2 ' ■ 



T a 



(C4) 



with new tensors r Q ° Q ° and (9 Q ° a °. 

We start by showing the gauge transformation ( C3 1 
In order to increase clarity we define 



3 



Qfc 

Sl — Sk 



Qai i^j^fc-ieifc 
[l]si ' "^[n]s t 

ai a,. 
P[l]si ' ' 'P[n]sk 



a k 

Sl — «* 



(C5) 



For this specific gauge proof the different Qui do not need 
to be of the same structure and neither do the pu\ . But 
in contrast to the p^ the still have to fulfill equation 
(§,i.e., 



— s> 



/3fc 

■Sl—Sh 

\a k _ 
-si---s k 



5°* , = «5' 9fcQfc , while 

- Sl-"Sfc 7^ ^(5l-"Sfe)(si-"Sft)5 



(C6) 



For fc — > oo this covers the first part of the heralded 
gauge equation ( C3 1 . The second part of equation ( C3 1 is 



proved by a straight forward application of the arguments 
used above on the right side of the MPS. 

Now, we have to show the MPS construction of the 
translational invariant ground state (C4|. As above, we 



assume that the ground state level of the Hamiltonian 
under consideration is g times degenerate due to broken 
translational invariance. Let Tl be the operator which 
shifts all sites of an MPS by one position to the left and 
Tr = {T L y l - For the MPS A representing any of the 
possible ground states we get 

T L -A £ A 
(T L ) 9 -A = A 

9-1 

t l -t = r, (en) 
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where T represents an unnormalized version of the trans- 
lational invariant ground state, which we like to con- 
struct. As an intermediate step we like to prove equation 
( C15 1 below. Therefor we have to look at the effect Tl 



has on A (using equation (CI I for A ) 



T I /n a -2°-i n a -i"i) \a a ^a Qi naia 2 
1 L [■ ■ ■W[ L ] S _ 1 '^[L]s ' A ' ^[R] S1 ' ^[R] s 2 



^Q_lCK_l ^)Q-iaO y^jCKQCKl 



[LI s _ 



'[R]s 



(C12) 



Next, we look at the MPS (T R ) 9 ■ T L ■ A, which has the 
form 



Q 



a fl _ 2Qg _ lAQg _ lQg _ 1( 



[L]s ' " * ^ [i] A ' ^[R]s fl "' 

(C13) 

Since the two MPS (Tr) s • T L • A = T L ■ A describe the 
same physical state we are allowed to apply the gauge 
transformation (CIO I and identify 



f[R] S„ s _ 



[i] so 



(C14) 

Inserting this expression into equation (C12 1 we arrive at 
the following description for Tl ■ A 



Appendix D: Davidson implementation 

In this subsection we introduce a practical implemen- 
tation resembling the Davidson [24] method based on re- 
cycled information of the previous round, which allows to 



improve the update equation (37) of the iterative eigen- 
vector solver presented in section |IV C| We adopt the 
same notation as in that section but mostly drop the in- 
dex [n] to keep the formulae clean. 

The best possible new vector |2t&+i) the iterative eigen- 
vector solver could come up with to replace equation ( 37 I 
is an orthonormalized version of |A ) = \E ) — |e ). 



fi.fleo) 
(Eo-l-i 

(e • I - 1 

(e • i — 5 



■\Eo) 
|A » 

■|A ) 
■|A ) 
■|A ) 
|A > 



E ■ \Eo) 

E ■ (|eo) + |Ao» 
- E • l) • |eo> 

-e -lV|e ) 



E 



(Dl) 



Tr 



Q 



a_ 2 a_ 1 
[L]*-i 



>a_ia 
{ [L] s 



•A 



'[R]S! 



-)QlQ!2 

Qa\a 2 
[R] s 2 ' 



(C15) 

As we see, applying the operator Tj, on A has the same 
effect as the replacement of the matrix A Q ° Q ° by w aoQ °. 
Since higher powers of Tl can be created by an iteration 
of the arguments just presented, it follows that the effect 
of any (Tl) 1 on A can be accounted for by an accordingly 
calculated ui"^ a °. Following that construction the only 

difference between the various MPS (Tl) 3 ■ A is their 



tensor Luf° a ° and hence the task of equation ( CI 1 1 to 



sum up these MPS reduces to a summation of the u)?S ao . 
In other words: Replacing the matrix A Q ° Q ° in the MPS 
A by T a ° a ° 



7=0 



(C16) 



results in the translational invariant MPS T. Of course, 
the same arguments can be applied to the MPS B giving 
us Tq = T = T q as claimed in equation ( C4 1 . 

Let us review our arguments: By virtue of equation 



(C16 I we can transform the MPS A and B given in equa- 



tion (CI) such that we end up with the translational 



invariant ground state Tq and T q as claimed in equation 
(C4|. Since we work under the condition that Tq = T q 



we are allowed to use the gauge transformation (CIO I to 



replace the q of Tq by the Q of Tq ■ The same replacement 
is possible in B because the q in B and T q are identical 
(as are the Q in A and Tq). This concludes the proof of 
equation (C2| we aimed for. 



where we used the definition (36 1 for \r). The Davidson 



method requires a workable approximation for the non 
trivial operator 2) = (Eq ■ I — Ilj . At this point we 

take advantage of the expectation that the operator 2)r n i 
calculated in round n should look pretty much the same 
as 2) r_ _i i calculated in round n — 1 



~ 25[n-l] 

-1 



n]0 



E 



n-l] 



(D2) 



Hence we use the accumulated data at the end of round 
ri—l for an efficient one time estimation of 5?[ n _i], which 
we will apply in round n. 

In order to calculate 2) we need a simplified form of 
H which allows easy inversion. We know k approximated 
eigenvectors |ej<fc) of H. In order to have an orthonormal 
basis for H we imagine N — k further |efc<j<jy), where 

dim ( H | — N x N. With that we approximate H as 



fc-i 



N 



^H|e l )(e i |+a-^|e J )( 

i—Q i—k 
fc-1 

2(1 -aj \e t )(e t \+a-I 



(D3) 



with an average eigenvalue a = const, for the unknown 
eigenvectors. One might be tempted to simplify equation 
(D3l using H|ej) ~ e i ■ |e») - but we do not, since the 



eigenvectors are not very well approximated except for 
[ep). To be able to perform the inversion in equation 



( D6 1 it suffices to resort to the exact result 
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Appendix E: Altered minimization 



— 6? • Si 



without 



(D4) 



which is a consequence of the construction (351. Further 



with the results gathered during the optimization ( IV C ) 
the H|ej) are as quickly calculated as the |e;). 



Next, we take the trace of equation (D3l and set a 
such that both sides are equal 



tr 



= (ei|H|e,) + a ■ (N - k) 



tr 



2-ii=0 e * 



N-k 



(D5) 



The trace of the exact H is efficiently calculated by al- 
ready tracing over its components L a ^ iiai and R a r^ a T 
before assembling them (151. 



Now we insert the approximated H (D3l in 53 
£ = (e ■ I - fif) 

« ((Eo-a)-I-J2(W-a) |e,:)( ei | ) .(D6) 



i=0 



The inversion is solved by 



fc-i 



a 



i=0 



E n - e,- 



(D7) 



as can be verified inserting the result in 53 ■ 53 _1 = I. 

The final question we have to answer is which value 
we assign to the unknown exact eigenvalue En. The best 
approximation (which we already used in equation ( D 1 1 ) 
is En ~ eo, but this produces a singularity in 53. There 
are two ways out. First, we can always pick En a little bit 
lower En := &a — s. Second, we should discard the trou- 
blesome term ~ |eo)(eo| in 53 anyway for the following 
reason: We replace equation (37 1 by 



where |r)_L|2li) 



|Bfc+i) = 53 ■ \r) 

refer l\ \J n ]\ 



Jn-1] 



that |eo)(eo|r) w • |eo) ■ Afterwards, residual parts of the 
|eo)(eo| term are exfiltrated again because |2l/c+i) has to 
be orthogonalized (and normalized) 



|2l 



fc+ly 



fc+ly 



(D9) 



These considerations are also part of the more elaborated 
Jacobi-Davidson [24] method to which this implementa- 
tion can be extended. 

We might further consider to omit terms with a 3> eo 
since their influence shrinks with (eo — e.;) _1 . At the end 
of the day, the effort to construct 53 as well as the effort 
for each application scale with N times the number of 
|ej)(ej| terms used in 53. 



Here we derive the missing equations of |IIID| First, 
we search for an approximation of \l/r] an d start by an 
alternative way of expressing them 



\ a P — n* Q T Ail 3 - \ a @ — A a i n*~iP 

A [L] — **[L] s ' J A [R] — A s ' W[R] s > 



(El) 



where we used the orthogonality Q of the Q and the 
decomposition (|8|. Next, we use the fact that the algo- 
rithm is tuned to produce consecutive tensors , 
which are quite similar. Hence we approximate the yet 
unknown Q^ L / R ^ (Ell by their known predecessor of the 



optimization round before, i.e. Q 



[n] * 
[L/R] 



Q 



[n—1] * 
[L/R] ■ 



l[n-l]* ay A[n\-yj3 
'[L]s A s 



T[n]a/3 _ A[n]a~f n [n-l]*J0 
A [R] — ' ^[R] s 



(E2) 



Next, we use the same idea and replace the 
definition (42 1 by the Q^/ R ^ 



10* 



with 



\*i)T 



X 



n] ck0 
[L]i 



<[R)s 



(E3) 



With that we like to calculate \A) = |2tj) - a, (44 1. There- 
for we first observe 



A 



n] a/3 



= A 



[L] 



(E4) 



(D8) 

) and with 1^) . ai 



where we used A — |2lj) • (41 1 in line two and equation 
(El ) in line three. Likewise we find X[R]i ■ a>i ^ A[#] and 



> I 21 *) + 4 (Q[L] ■ X [R] i + A [L] i ■ Q[R]) 

^\A) + -(Q lL] -X [ R ] + X lL] -Q [R] ) 

\\ A ) + \\ A ) + \Q\l] ■ \r] + ■ Q [R] 
\Q[l\ ■ A[l] + ^Q[r] ■ X [R] + 
+ ^Q[L] ■ A[r] + -X [L] ■ Q [R] 

2 (Q[L] ■ \syrn] + ^[sym] ' Q[R]) witn 
1 



Vsym] 



2 ( A [L] +\r]) i 



(E5) 



as claimed in equation ( 44 1 
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Appendix F: Mirror symmetry 



2. If 1. is fulfilled, then A? l ? r 



[rtl s 



(F2I 



Here we show that in case of a mirror symmetric 
Hamiltonian % i.e. a Hamiltonian that is invariant under 
inversion of the order of its sites 



U s l"' s 
tt Sl —s 



Hit 



(Fl) 



all tensors ^4j^" s r can be chosen mirror symmetrical in 
their auxiliary indices ai,a ri i.e. 



\n] s 



A \n\ s ' 



(F2) 



The very first H^j is constructed via the initialization 
procedure described in |IIE2| If we start with a mirror 
symmetric wave function and use Takagi's factorization 
as suggested in 3. of the initialization procedure, Hm is 
symmetric and the induction is well grounded. 

Proof for 1. All H represent a sum of operators ac- 



cording to equation (18 I 



This allows to impose an extra constraint on Ar n ] . Com- 
muting the indices a; , a r also results in 



\aia r 



[R] s 

\ OL r Ctl 

A [R] > 



(F3) 



as can be seen directly from the decomposition (|8]). Fur- 
ther it is possible to construct L a '^ iai and R<^a r 



such that they are identical. But therefor we need to re- 
sort to an alternative MPO construction for the Hamil- 
tonian, such that the MPO tensors of the left half are 
mirror symmetric to the tensors of the right half. This 
can be achieved if we include a special interface tensor in 
the middle where L a ^ iai and R a ,^ra r are connec ted to 
build HL This reduces the requirement in storage mem- 
ory roughly by a factor of 2, but has nearly no effect on 
the speed. Since storage capacities are usually not a big 
issue, we do not elaborate this point any further. 

One should be aware that equation (Fl ) enforces a mir- 



ror symmetric MPS. In case of a mirror symmetry break- 
ing ground state, the MPS will represent a superposition 
of both chiralitics, implying an unfavorably increased re- 
quirement in bond dimension. 

We further remark that our definition of a mirror sym- 
metric Hamiltonian does not forcedly imply a symmetry 
in real space. Although in practical application the order 
of the sites generally coincides with one specific spatial 
direction, there is no mathematical connection between 
the direction of space and the chosen order. 



Assuming a mirror symmetric Hamiltonian H (Fl I the 
claim of the mirror symmetric tensor AS" r (F2 1 can be 
proven iteratively: 



1. If A? 



A?r ai for all i 

h] s 



< n. then 



(181 is 



mirror symmetric, i.e. 



In] ss' 



In] ss' 



(F4) 



Each of these Hj contains an entire Hamiltonian MPO 
as central unit sandwiched by bra and ket MPS, with a 
hole where the new tensor ^4[ new ] = At n ] is supposed to 
be inserted. The Hamiltonian is guaranteed to be mirror 
symmetric, while there is no such condition for the MPS. 
The different MPS are encoded in the building blocks 
L"^' ' and R a rV r a r ^ which are constructed symmetri- 
cally (19 1, i.e., in contrast to the basic algorithm each 
new tensor is inserted in L aifJ - ia ' and R a r^- a r a t e q ua l 
footing. Hence, for each MPS exists a counterpart which 
contains exactly the same tensors in inverted order. In 
general, this alone is not enough, because mirror sym- 
metry also exchanges the left and right auxiliary indices. 
But, since all involved tensors ^4j*j" r = ^["' are sup- 
posed to be invariant under this kind of exchange (and 



Q^Ms = Q'lB^s (F3 1, as needed) , mirror symmetric coun- 



terparts for all involved MPS are guaranteed and with 
that EI is mirror symmetric. 

Proof for 2. The tensor -<4u l ! i" r is the result of the min- 
imization procedure described in more details in |IV C| 
Each element in this procedure maintains mirror sym- 



metry if . 



and the initial tensor A 



[refer] 



(27) are mir- 



ror symmetric. Since A^ fer ' is a superposition of mirror 
symmetric tensors, all conditions are met. 

Finally we remark that accumulating numerical errors 
might undermine the symmetry. Therefore we recom- 
mend to explicitly restore the symmetry of each -AS" r 
during its calculation. 



[1] U. Schollwoeck, Annals of Physics 326, 96 (2011). 

[2] S. R. White, Phys. Rev. Lett 69, 2863 (1992). 

[3] M. Fannes, B. Nachtergaele, and R. F. Werner, Comm. 

Math. Phys. 144, 443 (1992). 
[4] I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. 

Rev. Lett. 59, 799 (1987). 



[5] J. Dukelsky, M. A. Martin-delgado, T. Nishino, and 
G. Sierra, Europhys. Lett. 43, 457 (1998). 

[6] G. Vidal, Phys. Rev. Lett. 99, 220405 (2007). 

[7] R. Hiibener, C. Kruszynska, L. Hartmann, W. Diir, 
F. Verstraete, J. Eisert, and M. Plenio, Phys. Rev. A 
79, 022317 (2009). 



24 



[8] R. Hiibener, V. Nebendahl, and W. Diir, New J. Phys. 

12, 025004 (2010). [23 
[9] J. Eisert, Phys. Rev. Lett. 97, 260501 (2006). [24 
[10] N. Schuch, I. Cirac, and F. Verstraete, Phys. Rev. Lett. 

100, 250501 (2008). [25 
[11] F. J. Burnell, M. M. Parish, N. R. Cooper, and S. L. [26 

Sondhi, Phys. Rev. B 80, 174519 (2009). 
[12] J. Schachenmayer, I. Lesanovsky, A. Micheli, and A. J. [27 

Daley, New J. Phys. 12, 103044 (2010). 
[13] C. Trefzger, C. Menotti, B. Capogrosso-Sansone, and [28 

M. Lewenstein, J. Phys. B: At. Mol. Opt. Phys. 44, 

193001 (2011). [29 
[14] J. Cardy, Scaling and renormalization in statistical 

physics (Cambridge University Press, 1996). [30 
[15] B. Pirvu, G. Vidal, F. Verstraete, and L. Tagliacozzo, 

Phys. Rev. B 86, 075117 (2012). [31 
[16] A. Gendiar, M. Daniska, Y. Lee, and T. Nishino, Phys. 

Rev. A 83, 052118 (2011). [32 
[17] G. Vidal, Phys. Rev. Lett. 98, 070201 (2007). [33 
[18] I. P. McCulloch, J. Stat. Mech. , 10014 (2007). [34 
[19] V. Murg, J. Cirac, B. Pirvu, and F. Verstraete, New J. [35 

Phys. 12, 025012 (2010). 
[20] G. M. Crosswhite and D. Bacon, Phys. Rev. A 78, 012356 [36 

(2008). 

[21] F. Frowis, V. Nebendahl, and W. Diir, Phys. Rev. A 81, [37 
062337 (2010). 

[22] G. M. Crosswhite, A. C. Doherty, and G. Vidal, Phys. 



Rev. B 78, 035116 (2008). 

I. P. McCulloch, E-print: arXiv:0804.2509vl. 

G. L. G. Sleijpen and H. A. Van der Vorst, SIAM Review 
42, 267 (2000). 

M. B. Hastings, J STAT , 08024 (2007). 

H. Ueda, I. Maruyama, and K. Okunishi, J. Phys. Soc. 
Jpn. 80, 023001 (2011). 

H. Ueda and I. Maruyama, E-print: arXiv:1112.2075v2 
(2012). 

Y. Saad, Iterative methods for sparse linear systems ( 2nd 
ed.) (SIAM. ISBN 0898715342, 2003). 
C. Lanczos, Journal of research of the National Bureau 
of Standards 45, 255 (1951). 

W. Arnoldi, Quarterly of Applied Mathematics 9, 17 
(1951). 

S. Furukawa, M. Sato, S. Onoda, and A. Furusaki, Phys. 
Rev. B. 86, 094417 (2012). 

K. Okunishi, J. Phys. Soc. Jpn. 77, 114004 (2008). 
F. J. Burnell, E-print: arXiv:1004.5595v2. 
P. Bak and R. Bruinsma, Phys. Rev. Lett. 49, 249 (1982). 
T. Giamarchi, Quantum Physics in One Dimension (Ox- 
ford University Press, 2003). 

L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I. 
Latorre, Phys. Rev. B 78, 024410 (2008). 
F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. 
Moore, Phys. Rev. Lett. 102, 255701 (2009). 



